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Abstract 

We develop a theory for the nonequilibrium coherent transport through a 
mesoscopic region, based on the nonequilibrium Green function technique. 
The theory requires the weak coupling between the central mesoscopic region 
and the multiple electrodes connected to it, but allows arbitrary hopping and 
interaction in the central region. An equation determining the nonequilib- 
rium distribution in the central interacting region is derived and plays an 
important role in the theory. The theory is applied to two special cases for 
demonstrations, revealing the novel effects associated with the combination 
of phase coherence, Coulomb interaction, and nonequilibrium distribution. 
PACS numbers: 73.63.Kv, 85.35.Ds, 73.23.Hk, 73.40. Gk. 



In the realm of mesoscopic transport, the basic rules of the traditional electronics break- 
down. When the size of device is smaller than the electron mean free path in the material, 
electrons will behave like a wave rather than particles. Meanwhile, with the reduced size of 
device, the Coulomb interaction among electrons becomes important. An additional elec- 
tron has to overcome the repulsion from other electrons in the device before entering it. 
Moreover, due to lack of inelastic collisions, the thermal distribution in the device is no 
longer a equilibrium one when electrodes are biased with finite voltages. Therefore, a the- 
ory containing the above three "ingredients" — phase coherence, Coulomb interaction, and 
nonequilibrium — is of particular interest in the mesoscopic transport. 

Using nonequilibrium Keldysh formalism, Meir el al. derived a formula, in which the 
current flowing out of an electrode is expressed in term of the Green functions of the cen- 
tral region |Q. The remaining task is to find out these Green functions. Unfortunately, 
there seems no standard method to derive them, although various approximation schemes 
were used for special problem in specific system, e.g., large-N expansion ||, truncation for 
equation of motion ||, introducing an interpolative self-energy |4j], Ng's ansatz for lesser 
self-energy ||, etc. In fact, too much physics is hidden in the general Hamiltonian, and it is 
hopeless to invent a theory covering everything. If we restrict ourselves to the weak coupling 
case, the complex Kondo physics will be ruled out. In these circumstances, the Green func- 
tions in the "atomic limit" should be a good starting point for the construction of the full 
Green functions, and we shall show that a general solution to the problem is possible. We 
note that similar idea has been addressed in the linear response theory by several authors 
HI?], but hardly investigated in the nonlinear regime |§. 

The aim of this paper is to present a scheme for the calculation of the Green functions, 
and hence establish a theory of nonequilibrium coherent transport through an interacting 
mesoscopic region. As a price, the coupling between the central mesoscopic region and the 
electrodes is required to be relatively small. Electron transport through the mesoscopic 
region is viewed as a summation over various coherent processes via many-body quantum 
states, weighted by nonequilibrium thermal probabilities. The many-body quantum states 
in the central region can be found by exact diagonalization, while the nonequilibrium dis- 
tribution can be determined by an equation derived in the text . 

The mesoscopic system under consideration is modelled by the Hamiltonian (hereafter 
e = h = 1), H = H cent + Y,p H p + h t, where H cent = H cent ({ci}, {cj}) is a general 
Hamiltonian for the central region with M-sites (spin index has been absorbed into the 
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site index i), Hp = Ylp{ e k ~ Vp) a \k a Pk is the Hamiltonian for the /3th electrodes, and 
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Ht = Ylpki v Pi a 0k°i + H- c - * s the tunnel coupling between them. This Hamiltonian is ap- 



plicable to a large variety of mesoscopic systems, such as molecular devices, tunneling cou- 
pled carbon nanotubes, quantum dot arrays, Aharonov-Bohm rings embedded with quantum 
dots, etc. 

Define the Green function of the central region and the self-energy arise from the 
coupling with /5th electrodes as Glf' < (t 1 ,t 2 ) = ((ci(ti)|c](t 2 ))) r,a,< and ^f(t u t 2 ) = 

o' a ' vpj, where the superscript r, a, < denote for the retarded, ad- 
vanced, lesser Green function and self-energy, respectively; the subscript denotes for the 
Green function in the decoupling limit. Following Meir et al., the current flowing out of the 
/?th electrodes can be expressed in the compact form Ip = J ^Tr {2 Re [G(u;)S / 3(ci;)] < }, 
where [AB^ = A r B < + A < B a , G r,0,< (a;) and £ r,a,< (ci;) are the Fourier transformed Green 
function and self-energy. In the wide bandwidth limit, the self-energy can be evaluated 
as 5jg(o;) = — -I^ and S»(w) —ifp(uj)Tp, where Tp ti j = 2irDpVp i vpj with Dp being the 
density of states at the Fermi surface of the /3th electrode, fp(uS) = f{uj — Vp) with f(uj) 
being the Fermi distribution function. We assume that the coupling between the central 
region and the electrodes is relatively small, i.e., Tp^ <C ksT. Under this assumption, 
the Green function of the central region can approximately be written in a Dyson equation 
form G = g + gS G, with g = lim# T ^ G being the Green function in the decoupling limit, 
S = ^2p(£/3) being the self-energy contributed by the coupling with electrodes. (This ap- 
proximation is amount to a proper truncation of equation of motion.) Accordingly, the 
approximate "Dyson equation" for G r and the "Keldysh equation" for G< are 

Q r = r + rSS G r , (1) 

G< = G r S<G a . (2) 

Consequently, the current formula can be rewritten in a Landauer type as if in the nonin- 
teracting case M, 

J ^ = E/ ^-T a p(u)[fp(u)-f a (u)} , (3) 

with T a p(uj) = TrG r T a G a Tp. 

The remaining task is to calculate the retarded Green function in the decoupling limit. 
It is straightforward to exactly diagonalize H cent ({ci}, {c\}) in the particle occupation bases 
{( c m) Nm ■ • • ( c t) N2 ( c \) Nl 1°) where Ni = or 1, and obtain the 2 M eigenstates {E n ,\n)}. 
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Once again, under the weak coupling assumption, the density matrix operator of the central 
region is supposed to have the diagonal form p cent = J2 n Pn\ n ) (n\, with the constraint 
J2 n Pn = 1- Here, the central region is regarded as "system", while the electrodes are 
"environment" in local equilibrium. Given {P n }, the decoupled Green functions can be 
expressed as 

rem v ' 

=Y,2™PJ[w-(En-E m )]{m\A\n){n\B\m) , (5) 

nm 

where A and 5 are operators composed of {q} and {cj}. So the determination of the 
nonequilibrium distribution {P n } lies in the heart of the theory. In the linear response 
regime, i.e., \Vp — Vp\ <C ksT and hence Vp ~ Vb, the central region is in a thermal equi- 
librium, and the distribution can be written as P n = ^L e -{ E ^- N ^i)l k BT _ p or ^ cage Q f 
nonequilibrium, however, {P„} is determined by the coupling to electrodes with different 
chemical potentials, and in principle needs a self-consistent calculation. Our strategy is to 
choose a proper set of observables {0} and establish the equations of {P n } by the stationary 
condition 

(d t O) = (i[H,O]) = . (6) 

We point out that the 2 M conservables {Oi = \l) (|/) is the eigenstate of H cent ) are ideal 
candidates for the task. Notice that 

(i[H,O l }) = 2Re I ^EE [«[ci.O I ]| C }»^.«a /Jfc |at fc » t; /Ji ] < , (7) 

J /3k ij 

and make the approximation (([q, Oj]|cJ)) w (([q, 0;]|cj))o under the weak coupling approx- 
imation, one can derive a set of equations of {P n } as 

E E [PmMEn - E m ) - PME n - E m )] Yi m Q l nm = (I = 1, 2 • • • 2 M ) , (8) 

/3 nm 



where n and m run over all the eigenstates of H cent , T^ m = £V (m|cj|n) (n|c]|m)r /3i j : ,-, Q 



' nm 



— ^m«, and fp(uj) = 1 — fpiuj). Because J2i 10 (^1 = 1) the 2 M conservables can produce 
2 M -1 independent equations, and the constraint Y2 n Pn — 1 should be supplemented for 
completeness. Eq.(8) is the central result of this work, which determines the nonequilibrium 
distribution in an interacting system. With {P n } solved from the set of equations, one can 
calculate both nonequilibrium tunneling current and various quantities of the central region. 
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To sum up, Eq.(l), (3), (4), and (8) consist of the frame for the calculation of the 
nonequilibrium coherent transport through an interacting mesoscopic region, requiring weak 
coupling between the central region and the electrodes, but allowing arbitrary interaction 
and hopping in the central region. Below we shall apply the theory to two special cases for 
demonstrations. 

(1) a single quantum dot with multiple levels. Consider a single quantum dot (QD) 
connected to electrodes with finite bias voltages, which can be modelled by the Hamiltonian 

H cent = 22 EiHi + U niU i ' ( 9 ) 
i i<j 

where n« = cjcj is the particle number operator. For convenience, the particle occupation 
bases are numbered as F = Y^U Nf ■ T~ x and \F) = (c j M ) N M ■ ■ . (4)^(4)^ |0). Notice 
that H cent is already diagonalized in the {\F)} bases, and the eigenenergy of \F) is E\ F ) = 
YliEiNf + U^2 iK -N[ . The conservable \l) (/| can be written explicitly in the form of 
{cj} and {c\} as m l M ■ ■ ■m l 2 m l 1 where m\ = for N\ = 1 and m\ = 1 — for N\ = 0. The 
equations of {Pf} ar e simplified to 

M 

X)(-l)**r, [HE\ Fl) - E ]Fo) )P Fo - h(E lFl) - E ]Fo) )P Fl ] = (Z = 1, 2 ■ ■ ■ 2 M ) , (10) 

i=l 

where T { = ^ T Ph h^u) = J2p T^fpi"), ?m = !>,« = 2tt^|^| 2 , F 1 = l-N l i - + 2*" 1 , 
F = I — N\ ■ 2* _1 , and hi{uj) = 1 — hi{uj). The retarded Green function g r (a;) is obtained as 

r 1 

with Ef = Ei +U Ylj^i being the renormalized resonances. Specially, for M = 2, using 
(ni) = Pq\ + Pn, (n 2 ) = Pio + Pn, (^1^2) = Eq.(19) and Eq.(20) are equivalent to 

(n i ) = (l-(n 1 ))h i (E i ) + (n 1 )h i (E i + U) , (12) 
W = TT^rrhiE, + C/)(n 2 ) + -^-h 2 {E 2 + U)(n x ) , (13) 

1 1 -h 1 2 1 1 + J- 2 

with i = 3 — i for z = 1 or 2. Thus, our theory reproduces the correct results for the oc- 
cupation number (n») and the retarded Green function ((cj|c]))o in the limit of Tj — > [EJ, 
and derive the correlator {n\n 2 ) which is otherwise difficult to obtain. Fig.l shows the equi- 
librium and nonequilibrium distributions for the quantum dot containing two interacting 
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levels connected with two electrodes. One can see in the plot: (a) The complete Coulomb 
blockade in equilibrium is partially removed in nonequilibrium, i.e., the blockaded state can 
be occupied in some "windows" of the gate voltage, (b) The correlator (711712) is obviously 
unequal to (ni)(n2) in nonequilibrium, although approximately correct in equilibrium for 
nondegenerate levels, (c) The total occupation number has fractional steps in nonequilib- 
rium in contrast to integer steps in equilibrium, and the fluctuation of the total number is 
reminiscent of the shape of tunneling current. 

Next, we insert the QD to one arm of a Aharonov-Bohm (AB) interferometer f9j. The 
other arm (reference arm) is modelled by a quantum point contact, which can be described 
by adding the term ^2 k {We^a^ Lk aRk + H.c.) to H T jlOfl , with (f> being the AB phase induced 
by magnetic flux. Both the current formula and the equation of {P n } should be modified 
to include the "direct" coupling between electrodes, the details will be presented elsewhere. 
The background conductance of the reference arm (measured when QD is decoupled from 
both electrodes) is G = y (r+x) 2 w ^h x = n 2 W 2 D L D R . The effective conductance of QD is 
formally defined as Gdot = I{Vb, V g )/Vb — Gq. Fig. 2 shows the curves of Gdot in both linear 
and nonlinear transport regimes. Three features are noteworthy: (a) In linear transport, 
the conductance is contributed only by the ground state, while in nonlinear transport, the 
conductance is contributed by both ground and excited states, recognized by the substeps 
in the conductance plateaus and valleys, which is in agreement with the recent experiment 
pTl]] . (b) With the increase of the background conductance, Gdot vs V g curves exhibit Fano 
pattern, evolving from a peak (plateau) to a pair of peak and dip (plateau and valley), and 
finally to a dip (valley), which is originated from the constructive and destructive interference 
between a resonance and the uniform background, (c) The phase analysis shows that the 
dependence of Gdot on <fi across a resonance is quite different between linear and nonlinear 
transport. In linear transport, the dependence changes from cos0 to cos(0 + it) abruptly 
across a resonance, and cos(20) component only appears on the resonance. In nonlinear 
transport, cos(20) component is always accompanied with the cos0 component, and the 
crossover from cos0 to cos(</> + 71") occurs continuously. The "nonequilibrium- Fano" effect 
discussed here and the "Kondo-Fano" effect in [110] seem to share some similarities, although 
the mechanisms are basically different. 

(2) coupled double quantum dots. Consider two quantum dots coupled in series with left 
and right electrodes (N-QD=QD-N), each dot is capacitively coupled to a gate so that the 
energy level of the dot is tunable IT3]. The coupled double quantum dots can be modelled 
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by a 4-site Hamiltonian 

H cent = ^ ^ ^iffWiff + * X^ C !L C2,T + cL C ^) 

1=1,2 IT <7 

+ ^ U i n i} n il + ?7i 2 (wi T + ra u )(n 2T + n 2 +) . (15) 

i=l,2 

We first diagonalize H cent in the 2 4 particle occupation bases. Due to the particle number 
conservation and spin conservation, the 16 dimensional spaces can be divided into several 
subspaces 16 = 1 + (2 + 2) + (1 + 1 + 4) + (2 + 2) + 1, in which eigenstates are readily 
solved. In principle, one can find out 2 4 conservables written in {c ia } and {c\ a } as done in 
(1), although it is uneasy and unnecessary. We only need to calculate the effective coupling 
strength T^ m and put them into Eq.(8). With {P n } solved from the 2 4 linear equations, g r , 
hence G r and G< are available. Finally, the current formula in the N-QD=QD-N system 
can be simplified as 




We present the numerical results of the tunneling current I vs the resonant levels (Ei,E 2 ) 
in Fig.3. The bias voltage V b (V L = V b /2, V R = -V h /2) is fixed as 0.001, 0.3, and 0.6 for the 
graphs from top to bottom. The graph of V b = 0.001 is corresponding to the linear response 
regime, in which the thermal distribution in the central region is nearly equilibrium. Due 
to the intradot and interdot Coulomb interactions, the conductance peaks are arranged to a 
hexagon pattern, consistent with both experiment and the semiclassical Coulomb blockade 
model. The graph of V b = 0.6 is corresponding to the strong nonequilibrium case, in which 
thermal distribution of the central region is determined by two reservoirs with different 
chemical potentials. The effect of finite bias voltage on the tunneling current is two folded: 
(a) Pull each conductance peak along the direction of E\ — E 2 , and ultimately emerge them 
into a ridge. It is easy to understand the pulling effect by thinking of the noninteracting 
case (i.e., XJ\ — U 2 — U 12 = 0), where the conductance through N-QD=QD-N is allowed 
only when V L > E x « E 2 > V R . (b) Modify the hexagon pattern significantly, and break the 
symmetry with respect to E\ — E 2 . This can be attributed to the nonequilibrium occupation 
of the coupled quantum dots. To proceed, let us cut off the interdot hopping for a moment 
(i.e., t — 0). Then the electron filling to the left (right) dot is only related to the chemical 
potential of the left (right) electrode. The occupation configuration {(ni), (n 2 )} vs the 
energy levels (Ei,E 2 ) has the same shape of hexagon boundary as in the equilibrium case, 
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except for a displacement of (y , — ^). Therefore, the symmetry with respect to E\ = E 2 
is broken when V& ^ 0. Turning on the interdot hopping will make the problem much more 
complicated, either (ni) or (n 2 ) are not good quantum number, energy levels of each dot 
are hybridized into "molecular orbits" , and the occupation of the "molecule" is affected by 
both of the electrodes. Simple interpretation for this situation is beyond our intelligence. 

In conclusion, we have presented a theory dealing with nonequilibrium coherent transport 
through an arbitrary mesoscopic region, possibly containing strong Coulomb interactions. 
The only restriction of the theory is that the coupling between the central region and elec- 
trodes should be sufficient weak so that the central region may be regarded as a single 
quantum system. The key innovation of the theory is Eq. (8) which determines the nonequi- 
librium distribution in an interacting system. The general theory is applied to two special 
cases, a single quantum dot with multiple interacting levels and coupled double quantum 
dots, and novel behaviors are found in the nonlinear coherent transport. 

We would like to thank W. Li and Y. F. Yang for the valuable discussions. This project 
was supported by NSFC under Grant No. 10074001 and by the State Key Laboratory for 
Mesoscopic Physics in Peking University. 
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FIGURE CAPTIONS 



Fig. 1 Equilibrium and nonequilibrium distributions in a QD containing two energy levels 
indexed by % = 1,2, connected with two electrodes indexed by (3 = L,R. (a), (b) 
and (c) show the curves of (rii) (solid) and (n 2 ) (dotted), (nin 2 ) (solid) and {n\){n 2 ) 
(dotted), (n) = (m) + (n 2 ) (solid) and (Sn) = 4 [(n 2 ) - (n) 2 ] 1/2 (dotted) vs the gate 
voltage V g , respectively. The bias voltage H = Vl — V R is fixed at Vb = / Vb = 0.6 
for the left / right panel. Other parameters are: U — 1, k B T = 0.02, Ti L = Ti R = 
T 2 l = T 2R = 0.001, £?! = V g - 0.05, £ 2 = V g + 0.05. 

Fig. 2 Linear and nonlinear transport through an AB interferometer embedded with a 
QD (schematically shown in the inset), (a) and (b) show the curves of Gdot = 
I(Vb,V g )/Vb — G vs the gate voltage V g for = 0, with the bias voltage V& fixed 
at + and 0.6, respectively. The solid, dash, and dotted curves correspond to the 

2 2 

background conductance of the reference arm Go = 0, 0.5%- ; ir> respectively The QD 
contains three interacting levels, with level spacing AE = 0.1, interacting constant 
U = 1, and T L = T R = 0.001 < k B T = 0.02. (c) and (d) analyze the phase depen- 
dence of the G dot at the points marked in (a) and (b). Only the range of < < n is 
shown since G dot {<f)) = G dot (-(f)). 

Fig. 3 Surface graphs of tunneling current / vs the resonant levels (E\, E 2 ) in the N-QD=QD- 
N system, with sideview on the left and topview on the right. The bias voltage is fixed 
at V b = 0.001, 0.3, and 0.6 for the graphs from top to bottom. Other parameters 
are: U x = U 2 = 1, U 12 = 0.3, t = 0.1, k B T = 0.05, T L = T R = 0.01. (In the e- 
print version, Fig. 3 is separated into Fig. 3a, Fig. 3b and Fig. 3c corresponding to 
V b = 0.001, 0.3 and 0.6 respectively, to reduce the size of figures.) 
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